################################################################################
# PACKAGE AND DATA LOADING
################################################################################

pacman::p_load(fixest, readxl, modelsummary)

df_us <- readRDS('personnel_master_df.rds')

################################################################################
# SECTION 3: COEFFICIENT MAP
################################################################################

cm <- c(
  "wall_const"              = "Wall",
  "wall_const:mean_anxiety" = "Wall x Anxiety",
  "wall_const:mean_anger"   = "Wall x Anger",
  "mean_anxiety"            = "Anxiety",
  "mean_anger"              = "Anger",
  "staff"                   = "Border Staff",
  "staff_th"                = "Border Staff (thousands)",
  "log(staff)"              = "Border Staff (ln)",
  "southwest_border"        = "SW Border Staff",
  "southwest_border_th"     = "SW Border Staff (thousands)",
  "log(southwest_border)"   = "SW Border Staff (ln)",
  "bc_structure"            = "Border Crossing Infrastructure"
)

################################################################################
# SECTION 4: MODELS AND TABLES
################################################################################

# --- Table A.19: Personnel relationship with BIOs torture ---
m1 <- feols(log(amelia_fit+1) ~ staff_th,                                  data = df_us)
m2 <- feols(log(amelia_fit+1) ~ log(staff),                                data = df_us)
m3 <- feols(log(amelia_fit+1) ~ southwest_border_th,                       data = df_us)
m4 <- feols(log(amelia_fit+1) ~ wall_const + southwest_border_th,          data = df_us)
m5 <- feols(log(amelia_fit+1) ~ wall_const + log(southwest_border),        data = df_us)

models <- list("Torture (1)" = m1, "Torture (2)" = m2, "Torture (3)" = m3,
               "Torture (4)" = m4, "Torture (5)" = m5)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table A.19.docx',
             coef_map = cm, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors reported."))

# --- Table A.18: Wall and border infrastructure relationships with personnel ---
w1 <- feols(staff              ~ wall_const,   data = df_us)
w2 <- feols(log(staff)         ~ wall_const,   data = df_us)
w3 <- feols(southwest_border   ~ wall_const,   data = df_us)
w4 <- feols(log(southwest_border) ~ wall_const, data = df_us)
w5 <- feols(staff              ~ bc_structure, data = df_us)
w6 <- feols(log(staff)         ~ bc_structure, data = df_us)
w7 <- feols(southwest_border   ~ bc_structure, data = df_us)
w8 <- feols(log(southwest_border) ~ bc_structure, data = df_us)

models <- list("Staff"       = w1, "Staff (ln)"    = w2,
               "SW Staff"    = w3, "SW Staff (ln)" = w4,
               "Staff"       = w5, "Staff (ln)"    = w6,
               "SW Staff"    = w7, "SW Staff (ln)" = w8)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table A.18.docx',
             coef_map = cm, fmt = 1, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors reported."))